home *** CD-ROM | disk | FTP | other *** search
/ Amiga Plus 2000 #4 / Amiga Plus CD - 2000 - No. 4.iso / PowerPC / Dev / PPCRelease / Libm / vax / cabs.s < prev    next >
Encoding:
Text File  |  1993-09-23  |  4.5 KB  |  134 lines
  1. # Copyright (c) 1985 Regents of the University of California.
  2. # All rights reserved.
  3. #
  4. # Redistribution and use in source and binary forms, with or without
  5. # modification, are permitted provided that the following conditions
  6. # are met:
  7. # 1. Redistributions of source code must retain the above copyright
  8. #    notice, this list of conditions and the following disclaimer.
  9. # 2. Redistributions in binary form must reproduce the above copyright
  10. #    notice, this list of conditions and the following disclaimer in the
  11. #    documentation and/or other materials provided with the distribution.
  12. # 3. All advertising materials mentioning features or use of this software
  13. #    must display the following acknowledgement:
  14. #    This product includes software developed by the University of
  15. #    California, Berkeley and its contributors.
  16. # 4. Neither the name of the University nor the names of its contributors
  17. #    may be used to endorse or promote products derived from this software
  18. #    without specific prior written permission.
  19. #
  20. # THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  21. # ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  22. # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  23. # ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  24. # FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  25. # DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  26. # OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  27. # HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  28. # LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  29. # OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  30. # SUCH DAMAGE.
  31. #
  32. #    @(#)cabs.s    5.4 (Berkeley) 10/9/90
  33. #
  34.     .data
  35.     .align    2
  36. _sccsid:
  37. .asciz    "@(#)cabs.s    1.2 (Berkeley) 8/21/85; 5.4 (ucb.elefunt) 10/9/90"
  38.  
  39. # double precision complex absolute value
  40. # CABS by W. Kahan, 9/7/80.
  41. # Revised for reserved operands by E. LeBlanc, 8/18/82
  42. # argument for complex absolute value by reference, *4(ap)
  43. # argument for cabs and hypot (C fcns) by value, 4(ap)
  44. # output is in r0:r1 (error less than 0.86 ulps)
  45.  
  46.     .text
  47.     .align    1
  48.     .globl  _cabs
  49.     .globl  _hypot
  50.     .globl    _z_abs
  51.     .globl    libm$cdabs_r6
  52.     .globl    libm$dsqrt_r5
  53.  
  54. #    entry for c functions cabs and hypot
  55. _cabs:
  56. _hypot:
  57.     .word    0x807c        # save r2-r6, enable floating overflow
  58.     movq    4(ap),r0    # r0:1 = x
  59.     movq    12(ap),r2    # r2:3 = y
  60.     jmp    cabs2
  61. #    entry for Fortran use, call by:   d = abs(z)
  62. _z_abs:
  63.     .word    0x807c        # save r2-r6, enable floating overflow
  64.     movl    4(ap),r2    # indirect addressing is necessary here
  65.     movq    (r2)+,r0    # r0:1 = x
  66.     movq    (r2),r2        # r2:3 = y
  67.  
  68. cabs2:
  69.     bicw3    $0x7f,r0,r4    # r4 has signed biased exp of x
  70.     cmpw    $0x8000,r4
  71.     jeql    return        # x is a reserved operand, so return it
  72.     bicw3    $0x7f,r2,r5    # r5 has signed biased exp of y
  73.     cmpw    $0x8000,r5
  74.     jneq    cont        # y isn't a reserved operand
  75.     movq    r2,r0        # return y if it's reserved
  76.     ret
  77.  
  78. cont:
  79.     bsbb    regs_set    # r0:1 = dsqrt(x^2+y^2)/2^r6
  80.     addw2    r6,r0        # unscaled cdabs in r0:1
  81.     jvc    return        # unless it overflows
  82.     subw2    $0x80,r0    # halve r0 to get meaningful overflow
  83.     addd2    r0,r0        # overflow; r0 is half of true abs value
  84. return:
  85.     ret
  86.  
  87. libm$cdabs_r6:            # ENTRY POINT for cdsqrt
  88.                 # calculates a scaled (factor in r6)
  89.                 # complex absolute value
  90.  
  91.     movq    (r4)+,r0    # r0:r1 = x via indirect addressing
  92.     movq    (r4),r2        # r2:r3 = y via indirect addressing
  93.  
  94.     bicw3    $0x7f,r0,r5    # r5 has signed biased exp of x
  95.     cmpw    $0x8000,r5
  96.     jeql    cdreserved    # x is a reserved operand
  97.     bicw3    $0x7f,r2,r5    # r5 has signed biased exp of y
  98.     cmpw    $0x8000,r5
  99.     jneq    regs_set    # y isn't a reserved operand either?
  100.  
  101. cdreserved:
  102.     movl    *4(ap),r4    # r4 -> (u,v), if x or y is reserved
  103.     movq    r0,(r4)+    # copy u and v as is and return
  104.     movq    r2,(r4)        # (again addressing is indirect)
  105.     ret
  106.  
  107. regs_set:
  108.     bicw2    $0x8000,r0    # r0:r1 = dabs(x)
  109.     bicw2    $0x8000,r2    # r2:r3 = dabs(y)
  110.     cmpw    r0,r2
  111.     jgeq    ordered
  112.     movq    r0,r4
  113.     movq    r2,r0
  114.     movq    r4,r2        # force y's exp <= x's exp
  115. ordered:
  116.     bicw3    $0x7f,r0,r6    # r6 = exponent(x) + bias(129)
  117.     jeql    retsb        # if x = y = 0 then cdabs(x,y) = 0
  118.     subw2    $0x4780,r6    # r6 = exponent(x) - 14
  119.     subw2    r6,r0        # 2^14 <= scaled x < 2^15
  120.     bitw    $0xff80,r2
  121.     jeql    retsb        # if y = 0 return dabs(x)
  122.     subw2    r6,r2
  123.     cmpw    $0x3780,r2    # if scaled y < 2^-18
  124.     jgtr    retsb        #   return dabs(x)
  125.     emodd    r0,$0,r0,r4,r0    # r4 + r0:1 = scaled x^2
  126.     emodd    r2,$0,r2,r5,r2    # r5 + r2:3 = scaled y^2
  127.     addd2    r2,r0
  128.     addl2    r5,r4
  129.     cvtld    r4,r2
  130.     addd2    r2,r0        # r0:1 = scaled x^2 + y^2
  131.     jmp    libm$dsqrt_r5    # r0:1 = dsqrt(x^2+y^2)/2^r6
  132. retsb:
  133.     rsb            # error < 0.86 ulp
  134.